# Load Necessary Packages
library(glmnet)
library(foreign)
library(estimatr)
library(dplyr)
library(tidyr)
library(bcf)
library(dbarts)
library(texreg)
library(Hmisc)
library(grid)
library(gridExtra)
library(interplot)
library(margins)
library(car)
library(plotrix)
library(ggpubr)
library(effectsize)
library(MASS)
library(rms)
library(pwr)
library(DeclareDesign)
library(estimatr)
library(haven)
library(TOSTER)

# Load Data
setwd("/users/josephphillips/Dropbox/COVID-19/Data/")
data <- read.dta("US_pooled_complete_v12.dta")

# Load glmnet function
run_LASSO <- function(data, yvar){
  
  set.seed(2334)
  
  var_quo <- enquo(yvar)
  
  ## select x variables 
  dat <- data %>%
    dplyr::select(college, agegroup, male, married, frequent_church, region, pid3, ideo7, highincidence,crt, knowledge,nonwhite, polinterest,health_trust,media_trust)
  
  ## select y-variable
  out1 <- data %>% 
    dplyr::select(!!var_quo) %>% pull() # save as a vector instead of a data-frame to use model.matrix()
  
  merged1 <- cbind(out1, dat) %>% na.omit() ## merge and delete NA
  x1 <- model.matrix(out1 ~ ., data = merged1) ## make model matrix
  y1 <- merged1$out1 ## select y-variable from NA-deleted dataset
  mod1 <- cv.glmnet(x1, y1, alpha = 1, family = "gaussian") ## run lasso model
  coef <- coef(mod1) ## get coef  
  
  list(merged1, x1, y1, mod1, coef)
  
}


# Clean Controls
data$male <- as.numeric(data$male) -1
data$married <- as.numeric(data$married) - 1
data$frequent_church <- as.numeric(data$frequent_church) - 1
data$pid3 <- ifelse(data$democrat==1,0,ifelse(data$republican==1,2,ifelse(data$pid7=="Independent",1,NA)))
data$ideo7 <- as.numeric(data$ideo7)
data$high_incidence <- as.numeric(data$high_incidence) - 1
data$lives_highincidence <- as.numeric(data$lives_highincidence) - 1
data$nonwhite <- as.numeric(data$nonwhite) - 1
data$polinterest <- as.numeric(data$polinterest)
data$factcheck_treatment <- as.numeric(data$factcheck_treatment) - 1
data$W3factcheck_treatment <- as.numeric(data$W3factcheck_treatment) -1
data$approve_trmp <- as.numeric(data$approve_trmp) - 1
data$conspiracy <- data$conspiracy-1
data$health_trust <- data$health_trust-1
data$media_trust <- data$media_trust-1
data$therm_trump <- data$therm_trump/100

# Clean DVs
data$wear_mask <- 6-(data$wear_mask) ## It was reverse-coded
data$mask_effective <- as.numeric(data$mask_effective)
## Note: with affective polarization, the party affect variables are actually mislabeled
data$affective_polarization <- ifelse(data$pid3==0,data$feeling_Rep-data$feeling_Dem,ifelse(data$pid3==2,data$feeling_Dem-data$feeling_Rep,NA))

# Clean Treatment
data$american_treatment <- ifelse(data$exposure_treat=="Americans",1,0)
data$copartisan_treatment <- ifelse((data$pid3==0 & data$exposure_treat=="Democrats") | (data$pid3==2 & data$exposure_treat=="Republicans"),1,0)
data$contrapartisan_treatment <- ifelse((data$pid3==0 & data$exposure_treat=="Republicans") | (data$pid3==2 & data$exposure_treat=="Democrats"),1,0)

# Create underestimation/overestimation in mask wearing
data$mask_allthetime <- ifelse(data$W2_covid19_actions_2=="Most of the time" | data$W2_covid19_actions_2=="All of the time",1,0)
data$underestimates_masks_american <- ifelse(data$wear_mask_pub<70,1,0)
data$overestimates_masks_american <- ifelse(data$wear_mask_pub>=90,1,0)
data$underestimates_masks_copartisans <- ifelse(data$pid3==2 & data$wear_mask_Rep<47,1,ifelse(data$pid3==0 & data$wear_mask_Dem<80,1,ifelse(data$pid3==1,NA,0)))
data$overestimates_masks_copartisans <- ifelse(data$pid3==2 & data$wear_mask_Rep>56,1,ifelse(data$pid3==0 & data$wear_mask_Dem>98,1,ifelse(data$pid3==1,NA,0)))
data$underestimates_masks_contrapartisans <- ifelse(data$pid3==0 & data$wear_mask_Rep<47,1,ifelse(data$pid3==2 & data$wear_mask_Dem<80,1,ifelse(data$pid3==1,NA,0)))
data$overestimates_masks_contrapartisans <- ifelse(data$pid3==0 & data$wear_mask_Rep>56,1,ifelse(data$pid3==2 & data$wear_mask_Dem>98,1,ifelse(data$pid3==1,NA,0))) ## Note: It's impossible to overestimate Democrat mask wearing by 10%

data$underestimates_masks_republicans <- ifelse(data$wear_mask_Rep<47,1,0)
data$overestimates_masks_republicans <- ifelse(data$wear_mask_Rep>56,1,0)
data$underestimates_masks_democrats <- ifelse(data$wear_mask_Dem<80,1,0)
data$overestimates_masks_democrats <- ifelse(data$wear_mask_Dem>98,1,0)

data$treatment <- as.factor(ifelse(data$american_treatment==1,"American",ifelse(data$copartisan_treatment==1,"Co-Partisan",ifelse(data$contrapartisan_treatment==1,"Out-Partisan","Control"))))
data.partisans <- subset(data,pid3==0 | pid3==2)
data.full <- subset(data,treatment=="American" | treatment=="Control")

# Power Analysis
## Model 1: Baseline Model, Mask Wearing
mod1.1 <- summary(lm_robust(wear_mask~american_treatment+copartisan_treatment+contrapartisan_treatment+male+pid3+ideo7+health_trust,data=data.partisans))
mod1.2 <- summary(lm_robust(wear_mask~american_treatment+male+pid3+ideo7+health_trust,data=data.full))
## Model 2: Underestimation Model, Mask Wearing
mod2.1 <- summary(lm_robust(wear_mask~(american_treatment*underestimates_masks_american)+(copartisan_treatment*underestimates_masks_copartisans)+(contrapartisan_treatment*underestimates_masks_contrapartisans)+(american_treatment*overestimates_masks_american)+(copartisan_treatment*overestimates_masks_copartisans)+(contrapartisan_treatment*overestimates_masks_contrapartisans)+male+pid3+ideo7+health_trust,data=data.partisans))
mod2.2 <- summary(lm_robust(wear_mask~(american_treatment*underestimates_masks_american)+(american_treatment*overestimates_masks_american)+male+pid3+ideo7+health_trust,data=data.full))
## Model 3: Baseline Model, Mask Effectiveness
mod3.1 <- summary(lm_robust(mask_effective~american_treatment+copartisan_treatment+contrapartisan_treatment+pid3+ideo7+health_trust+media_trust,data=data.partisans))
mod3.2 <- summary(lm_robust(mask_effective~american_treatment+pid3+ideo7+health_trust+media_trust,data=data.full))
## Model 4: Partisan Model, Mask Wearing
mod4 <- summary(lm_robust(wear_mask~(american_treatment*republican)+(copartisan_treatment*republican)+(contrapartisan_treatment*republican)+male+ideo7+health_trust,data=data.partisans))
## Model 5: Partisan Model, Mask Effectiveness
mod5 <- summary(lm_robust(mask_effective~(american_treatment*republican)+(copartisan_treatment*republican)+(contrapartisan_treatment*republican)+ideo7+health_trust+media_trust,data=data.partisans))
## Model 6: Fact-Check Model, Mask Wearing
mod6 <- summary(lm_robust(wear_mask~(american_treatment*factcheck_treatment)+(american_treatment*W3factcheck_treatment)+(copartisan_treatment*factcheck_treatment)+(copartisan_treatment*W3factcheck_treatment)+(contrapartisan_treatment*factcheck_treatment)+(contrapartisan_treatment*W3factcheck_treatment)+male+pid3+ideo7+health_trust,data=data.partisans))
## Model 7: Fact-Check Model, Mask Effectiveness
mod7 <- summary(lm_robust(mask_effective~(american_treatment*factcheck_treatment)+(american_treatment*W3factcheck_treatment)+(copartisan_treatment*factcheck_treatment)+(copartisan_treatment*W3factcheck_treatment)+(contrapartisan_treatment*factcheck_treatment)+(contrapartisan_treatment*W3factcheck_treatment)+pid3+ideo7+health_trust+media_trust,data=data.partisans))
## Model 8: Baseline Model, Affective Polarization
mod8 <- summary(lm_robust(affective_polarization~american_treatment+copartisan_treatment+contrapartisan_treatment+ideo7+polinterest,data=data.partisans))
## Model 9: Partisan Model, Affective Polarization
mod9 <- summary(lm_robust(affective_polarization~(american_treatment*republican)+(copartisan_treatment*republican)+(contrapartisan_treatment*republican)+ideo7+polinterest,data=data.partisans))

## Set up Power Analysis N's
n.partisans <- 2519
n.full <- 1709

## SD of residuals for all models
sd.1.1 <- sqrt(mod1.1$res_var)
sd.1.2 <- sqrt(mod1.2$res_var)
sd.2.1 <- sqrt(mod2.1$res_var)
sd.2.2 <- sqrt(mod2.2$res_var)
sd.3.1 <- sqrt(mod3.1$res_var)
sd.3.2 <- sqrt(mod3.2$res_var)
sd.4 <- sqrt(mod4$res_var)
sd.5 <- sqrt(mod5$res_var)
sd.6 <- sqrt(mod6$res_var)
sd.7 <- sqrt(mod7$res_var)
sd.8 <- sqrt(mod8$res_var)
sd.9 <- sqrt(mod9$res_var)

## Model ATE's
### Baseline Mask-Wearing (Partisans)
ate_a_1.1 <- mod1.1$coefficients[2]
ate_b_1.1 <- mod1.1$coefficients[3]
ate_c_1.1 <- mod1.1$coefficients[4]
### Baseline Mask-Wearing (Full)
ate_a_1.2 <- mod1.2$coefficients[2]
### Underestimation Mask-Wearing (Partisans)
ate_a_2.1 <- mod2.1$coefficients[2]
ate_b_2.1 <- mod2.1$coefficients[4]
ate_c_2.1 <- mod2.1$coefficients[6]
cov_1_2.1 <- mod2.1$coefficients[3]
cov_3_2.1 <- mod2.1$coefficients[5]
cov_5_2.1 <- mod2.1$coefficients[7]
cov_2_2.1 <- mod2.1$coefficients[8]
cov_4_2.1 <- mod2.1$coefficients[9]
cov_6_2.1 <- mod2.1$coefficients[10]
ate_a_cov1_2.1<- mod2.1$coefficients[15]
ate_a_cov2_2.1<- mod2.1$coefficients[18]
ate_b_cov3_2.1<- mod2.1$coefficients[16]
ate_b_cov4_2.1<- mod2.1$coefficients[19]
ate_c_cov5_2.1<- mod2.1$coefficients[17]
ate_c_cov6_2.1<- mod2.1$coefficients[20]
### Underestimation Mask-Wearing (Full)
ate_a_2.2 <- mod2.2$coefficients[2]
cov_1_2.2 <- mod2.2$coefficients[3]
cov_2_2.2 <- mod2.2$coefficients[4]
ate_a_cov1_2.2 <- mod2.2$coefficients[9]
ate_a_cov2_2.2 <- mod2.2$coefficients[10]
### Baseline Mask Effectiveness (Partisans)
ate_a_3.1 <- mod3.1$coefficients[2]
ate_b_3.1 <- mod3.1$coefficients[3]
ate_c_3.1 <- mod3.1$coefficients[4]
### Baseline Mask Effectiveness (Full)
ate_a_3.2 <- mod3.2$coefficients[2]
### Partisan Mask-Wearing
ate_a_4 <- mod4$coefficients[2]
ate_b_4 <- mod4$coefficients[4]
ate_c_4 <- mod4$coefficients[5]
cov_1_4 <- mod4$coefficients[3]
ate_a_cov1_4 <- mod4$coefficients[9]
ate_b_cov1_4 <- mod4$coefficients[10]
ate_c_cov1_4 <- mod4$coefficients[11]
### Partisan Mask-Effectiveness
ate_a_5 <- mod5$coefficients[2]
ate_b_5 <- mod5$coefficients[4]
ate_c_5 <- mod5$coefficients[5]
cov_1_5 <- mod5$coefficients[3]
ate_a_cov1_5 <- mod5$coefficients[9]
ate_b_cov1_5 <- mod5$coefficients[10]
ate_c_cov1_5 <- mod5$coefficients[11]
### Fact-Check Mask-Wearing
ate_a_6 <- mod6$coefficients[2]
ate_b_6 <- mod6$coefficients[5]
ate_c_6 <- mod6$coefficients[6]
cov_1_6 <- mod6$coefficients[3]
cov_2_6 <- mod6$coefficients[4]
ate_a_cov1_6 <- mod6$coefficients[11]
ate_a_cov2_6 <- mod6$coefficients[12]
ate_b_cov1_6 <- mod6$coefficients[13]
ate_b_cov2_6 <- mod6$coefficients[14]
ate_c_cov1_6 <- mod6$coefficients[15]
ate_c_cov2_6 <- mod6$coefficients[16]
### Fact-Check Mask Effectiveness
ate_a_7 <- mod7$coefficients[2]
ate_b_7 <- mod7$coefficients[5]
ate_c_7 <- mod7$coefficients[6]
cov_1_7 <- mod7$coefficients[3]
cov_2_7 <- mod7$coefficients[4]
ate_a_cov1_7 <- mod7$coefficients[11]
ate_a_cov2_7 <- mod7$coefficients[12]
ate_b_cov1_7 <- mod7$coefficients[13]
ate_b_cov2_7 <- mod7$coefficients[14]
ate_c_cov1_7 <- mod7$coefficients[15]
ate_c_cov2_7 <- mod7$coefficients[16]
### Baseline Affective Polarization
ate_a_8 <- mod8$coefficients[2]
ate_b_8 <- mod8$coefficients[3]
ate_c_8 <- mod8$coefficients[4]
### Partisan Affective Polarization
ate_a_9 <- mod9$coefficients[2]
ate_b_9 <- mod9$coefficients[4]
ate_c_9 <- mod9$coefficients[5]
cov_1_9 <- mod9$coefficients[3]
ate_a_cov1_9 <- mod9$coefficients[8]
ate_b_cov1_9 <- mod9$coefficients[9]
ate_c_cov1_9 <- mod9$coefficients[10]

## Declare Populations
set.seed(1000)
population1.1 <- declare_population(N=n.partisans,ate_A=ate_a_1.1,ate_B=ate_b_1.1,ate_C=ate_c_1.1,
                                    u=rnorm(n.partisans,mean=0,sd=sd.1.1))
set.seed(1000)
population1.2 <- declare_population(N=n.full,ate_A=ate_a_1.2,
                                    u=rnorm(n.full,mean=0,sd=sd.1.2))
set.seed(1000)
population2.1 <- declare_population(N=n.partisans,ate_A=ate_a_2.1,ate_B=ate_b_2.1,ate_C=ate_c_2.1,covariate1=cov_1_2.1,covariate2=cov_2_2.1,covariate3=cov_3_2.1,covariate4=cov_4_2.1,covariate5=cov_5_2.1,covariate6=cov_6_2.1,
                                    ate_A_covariate1=ate_a_cov1_2.1,ate_A_covariate2=ate_a_cov2_2.1,ate_B_covariate3=ate_b_cov3_2.1,ate_B_covariate4=ate_b_cov4_2.1,ate_C_covariate5=ate_c_cov5_2.1,ate_C_covariate6=ate_c_cov6_2.1,
                                    u=rnorm(n.partisans,mean=0,sd=sd.2.1),
                                    observed_cov1=data.partisans$underestimates_masks_american,
                                    observed_cov2=data.partisans$overestimates_masks_american,
                                    observed_cov3=data.partisans$underestimates_masks_copartisans,
                                    observed_cov4=data.partisans$overestimates_masks_copartisans,
                                    observed_cov5=data.partisans$underestimates_masks_contrapartisans,
                                    observed_cov6=data.partisans$overestimates_masks_contrapartisans)
set.seed(1000)
population2.2 <- declare_population(N=n.full,ate_A=ate_a_2.2,covariate1=cov_1_2.2,covariate2=cov_2_2.2,
                                    ate_A_covariate1=ate_a_cov1_2.2,ate_A_covariate2=ate_a_cov2_2.2,
                                    u=rnorm(n.full,mean=0,sd=sd.2.2),
                                    observed_cov1=data.full$underestimates_masks_american,
                                    observed_cov2=data.full$overestimates_masks_american)
set.seed(1000)
population3.1 <- declare_population(N=n.partisans,ate_A=ate_a_3.1,ate_B=ate_b_3.1,ate_C=ate_c_3.1,
                                    u=rnorm(n.partisans,mean=0,sd=sd.3.1))
set.seed(1000)
population3.2 <- declare_population(N=n.full,ate_A=ate_a_3.2,
                                    u=rnorm(n.full,mean=0,sd=sd.3.2))
set.seed(1000)
population4 <- declare_population(N=n.partisans,ate_A=ate_a_4,ate_B=ate_b_4,ate_C=ate_c_4,covariate1=cov_1_4,
                                  ate_A_covariate1=ate_a_cov1_4,ate_B_covariate1=ate_b_cov1_4,ate_C_covariate1=ate_c_cov1_4,
                                  u=rnorm(n.partisans,mean=0,sd=sd.4),
                                  observed_cov1=data.partisans$republican)
set.seed(1000)
population5 <- declare_population(N=n.partisans,ate_A=ate_a_5,ate_B=ate_b_5,ate_C=ate_c_5,covariate1=cov_1_5,
                                  ate_A_covariate1=ate_a_cov1_5,ate_B_covariate1=ate_b_cov1_5,ate_C_covariate1=ate_c_cov1_5,
                                  u=rnorm(n.partisans,mean=0,sd=sd.5),
                                  observed_cov1=data.partisans$republican)
set.seed(1000)
population6 <- declare_population(N=n.partisans,ate_A=ate_a_6,ate_B=ate_b_6,ate_C=ate_c_6,covariate1=cov_1_6,covariate2=cov_2_6,
                                  ate_A_covariate1=ate_a_cov1_6,ate_A_covariate2=ate_a_cov2_6,ate_B_covariate1=ate_b_cov1_6,ate_B_covariate2=ate_b_cov2_6,ate_C_covariate1=ate_c_cov1_6,ate_C_covariate2=ate_c_cov2_6,
                                  u=rnorm(n.partisans,mean=0,sd=sd.6),
                                  observed_cov1=data.partisans$factcheck_treatment,
                                  observed_cov2=data.partisans$W3factcheck_treatment)
set.seed(1000)
population7 <- declare_population(N=n.partisans,ate_A=ate_a_7,ate_B=ate_b_7,ate_C=ate_c_7,covariate1=cov_1_7,covariate2=cov_2_7,
                                  ate_A_covariate1=ate_a_cov1_7,ate_A_covariate2=ate_a_cov2_7,ate_B_covariate1=ate_b_cov1_7,ate_B_covariate2=ate_b_cov2_7,ate_C_covariate1=ate_c_cov1_7,ate_C_covariate2=ate_c_cov2_7,
                                  u=rnorm(n.partisans,mean=0,sd=sd.7),
                                  observed_cov1=data.partisans$factcheck_treatment,
                                  observed_cov2=data.partisans$W3factcheck_treatment)
set.seed(1000)
population8 <- declare_population(N=n.partisans,ate_A=ate_a_8,ate_B=ate_b_8,ate_C=ate_c_8,
                                  u=rnorm(n.partisans,mean=0,sd=sd.8))
set.seed(1000)
population9 <- declare_population(N=n.partisans,ate_A=ate_a_9,ate_B=ate_b_9,ate_C=ate_c_9,covariate1=cov_1_9,
                                  ate_A_covariate1=ate_a_cov1_9,ate_B_covariate1=ate_b_cov1_9,ate_C_covariate1=ate_c_cov1_9,
                                  u=rnorm(n.partisans,mean=0,sd=sd.9),
                                  observed_cov1=data.partisans$republican)

## Potential Outcomes
potential_1.1 <- declare_potential_outcomes(Y_Z_0=u,
                                            Y_Z_1=u+ate_A,
                                            Y_Z_2=u+ate_B,
                                            Y_Z_3=u+ate_C)

potential_1.2 <- declare_potential_outcomes(Y_Z_0=u,
                                            Y_Z_1=u+ate_A)

potential_2.1 <- declare_potential_outcomes(Y_Z_0=u+covariate1*observed_cov1+covariate2*observed_cov2+covariate3*observed_cov3+covariate4*observed_cov4+covariate5*observed_cov5+covariate6*observed_cov6,
                                            Y_Z_1=u+covariate1*observed_cov1+covariate2*observed_cov2+ate_A+ate_A_covariate1*observed_cov1+ate_A_covariate2*observed_cov2,
                                            Y_Z_2=u+covariate3*observed_cov3+covariate4*observed_cov4+ate_B+ate_B_covariate3*observed_cov3+ate_B_covariate4*observed_cov4,
                                            Y_Z_3=u+covariate5*observed_cov5+covariate6*observed_cov6+ate_C+ate_C_covariate5*observed_cov5+ate_C_covariate6*observed_cov6)

potential_2.2 <- declare_potential_outcomes(Y_Z_0=u+covariate1*observed_cov1+covariate2*observed_cov2,
                                            Y_Z_1=u+covariate1*observed_cov1+covariate2*observed_cov2+ate_A+ate_A_covariate1*observed_cov1+ate_A_covariate2*observed_cov2)

potential_3.1 <- declare_potential_outcomes(Y_Z_0=u,
                                            Y_Z_1=u+ate_A,
                                            Y_Z_2=u+ate_B,
                                            Y_Z_3=u+ate_C)

potential_3.2 <- declare_potential_outcomes(Y_Z_0=u,
                                            Y_Z_1=u+ate_A)

potential_4 <- declare_potential_outcomes(Y_Z_0=u+covariate1*observed_cov1,
                                          Y_Z_1=u+covariate1*observed_cov1+ate_A+ate_A_covariate1*observed_cov1,
                                          Y_Z_2=u+covariate1*observed_cov1+ate_B+ate_B_covariate1*observed_cov1,
                                          Y_Z_3=u+covariate1*observed_cov1+ate_C+ate_C_covariate1*observed_cov1)

potential_5  <- declare_potential_outcomes(Y_Z_0=u+covariate1*observed_cov1,
                                           Y_Z_1=u+covariate1*observed_cov1+ate_A+ate_A_covariate1*observed_cov1,
                                           Y_Z_2=u+covariate1*observed_cov1+ate_B+ate_B_covariate1*observed_cov1,
                                           Y_Z_3=u+covariate1*observed_cov1+ate_C+ate_C_covariate1*observed_cov1)

potential_6 <- declare_potential_outcomes(Y_Z_0=u+covariate1*observed_cov1+covariate2*observed_cov2,
                                          Y_Z_1=u+covariate1*observed_cov1+covariate2*observed_cov2+ate_A+ate_A_covariate1*observed_cov1+ate_A_covariate2*observed_cov2,
                                          Y_Z_2=u+covariate1*observed_cov1+covariate2*observed_cov2+ate_B+ate_B_covariate1*observed_cov1+ate_B_covariate2*observed_cov2,
                                          Y_Z_3=u+covariate1*observed_cov1+covariate2*observed_cov2+ate_C+ate_C_covariate1*observed_cov1+ate_C_covariate2*observed_cov2)

potential_7 <- declare_potential_outcomes(Y_Z_0=u+covariate1*observed_cov1+covariate2*observed_cov2,
                                          Y_Z_1=u+covariate1*observed_cov1+covariate2*observed_cov2+ate_A+ate_A_covariate1*observed_cov1+ate_A_covariate2*observed_cov2,
                                          Y_Z_2=u+covariate1*observed_cov1+covariate2*observed_cov2+ate_B+ate_B_covariate1*observed_cov1+ate_B_covariate2*observed_cov2,
                                          Y_Z_3=u+covariate1*observed_cov1+covariate2*observed_cov2+ate_C+ate_C_covariate1*observed_cov1+ate_C_covariate2*observed_cov2)

potential_8 <- declare_potential_outcomes(Y_Z_0=u,
                                          Y_Z_1=u+ate_A,
                                          Y_Z_2=u+ate_B,
                                          Y_Z_3=u+ate_C)

potential_9 <- declare_potential_outcomes(Y_Z_0=u+covariate1*observed_cov1,
                                          Y_Z_1=u+covariate1*observed_cov1+ate_A+ate_A_covariate1*observed_cov1,
                                          Y_Z_2=u+covariate1*observed_cov1+ate_B+ate_B_covariate1*observed_cov1,
                                          Y_Z_3=u+covariate1*observed_cov1+ate_C+ate_C_covariate1*observed_cov1)

## Create Subpopulations
underestimates_americans_2.1 <- population2.1() %>% filter(observed_cov1==1)
overestimates_americans_2.1 <- population2.1() %>% filter(observed_cov2==1)
underestimates_copartisans_2.1 <- population2.1() %>% filter(observed_cov3==1)
overestimates_copartisans_2.1 <- population2.1() %>% filter(observed_cov4==1)
underestimates_contrapartisans_2.1 <- population2.1() %>% filter(observed_cov5==1)
overestimates_contrapartisans_2.1 <- population2.1() %>% filter
underestimates_americans_2.2 <- population2.2() %>% filter(observed_cov1==1)
overestimates_americans_2.2 <- population2.2() %>% filter(observed_cov2==1)
republican_4 <- population4() %>% filter(observed_cov1==1)
republican_5 <- population5() %>% filter(observed_cov1==1)
factcheck_6 <- population6() %>% filter(observed_cov1==1)
W3factcheck_6 <- population6() %>% filter(observed_cov2==1)
factcheck_7 <- population7() %>% filter(observed_cov1==1)
W3factcheck_7 <- population7() %>% filter(observed_cov2==1)
republican_9 <- population9() %>% filter(observed_cov1==1)

## Randomizations
randomization_partisans <- declare_assignment(conditions=0:3,legacy=T)
randomization_full <- declare_assignment(conditions=0:1,legacy=T)

## Define Estimands
estimands_1.1 <- declare_inquiry(american=mean(Y_Z_1-Y_Z_0-ate_A),
                                 copartisans=mean(Y_Z_2-Y_Z_0-ate_B),
                                 outpartisans=mean(Y_Z_3-Y_Z_0-ate_C))

estimands_1.2 <- declare_inquiry(american=mean(Y_Z_1-Y_Z_0-ate_A))

estimands_2.1.accurate_a <- declare_inquiry(american_accuratenat=mean(Y_Z_1-Y_Z_0-ate_A),
                                            copartisans_accuratenat=mean(Y_Z_2-Y_Z_0-ate_B),
                                            outpartisans_accuratenat=mean(Y_Z_3-Y_Z_0-ate_C),
                                            subset=(observed_cov1==0 & observed_cov2==0))

estimands_2.1.underestimate_a <- declare_inquiry(american_underestimatenat=mean(Y_Z_1-Y_Z_0-ate_A),
                                                 copartisans_underestimatenat=mean(Y_Z_2-Y_Z_0-ate_B),
                                                 outpartisans_underestimatenat=mean(Y_Z_3-Y_Z_0-ate_C),
                                                 subset=(observed_cov1==1))

estimands_2.1.overestimate_a <- declare_inquiry(american_overestimatenat=mean(Y_Z_1-Y_Z_0-ate_A),
                                                copartisans_overestimatenat=mean(Y_Z_2-Y_Z_0-ate_B),
                                                outpartisans_overestimatenat=mean(Y_Z_3-Y_Z_0-ate_C),
                                                subset=(observed_cov2==1))

estimands_2.1.accurate_c <- declare_inquiry(american_accurateco=mean(Y_Z_1-Y_Z_0-ate_A),
                                            copartisans_accurateco=mean(Y_Z_2-Y_Z_0-ate_B),
                                            outpartisans_accurateco=mean(Y_Z_3-Y_Z_0-ate_C),
                                            subset=(observed_cov3==0 & observed_cov4==0))

estimands_2.1.underestimate_c <- declare_inquiry(american_underestimateco=mean(Y_Z_1-Y_Z_0-ate_A),
                                                 copartisans_underestimateco=mean(Y_Z_2-Y_Z_0-ate_B),
                                                 outpartisans_underestimateco=mean(Y_Z_3-Y_Z_0-ate_C),
                                                 subset=(observed_cov3==1))

estimands_2.1.overestimate_c <- declare_inquiry(american_overestimateco=mean(Y_Z_1-Y_Z_0-ate_A),
                                                copartisans_overestimateco=mean(Y_Z_2-Y_Z_0-ate_B),
                                                outpartisans_overestimateco=mean(Y_Z_3-Y_Z_0-ate_C),
                                                subset=(observed_cov4==1))

estimands_2.1.accurate_o <- declare_inquiry(american_accurateout=mean(Y_Z_1-Y_Z_0-ate_A),
                                            copartisans_accurateout=mean(Y_Z_2-Y_Z_0-ate_B),
                                            outpartisans_accurateout=mean(Y_Z_3-Y_Z_0-ate_C),
                                            subset=(observed_cov5==0 & observed_cov6==0))

estimands_2.1.underestimate_o <- declare_inquiry(american_underestimateout=mean(Y_Z_1-Y_Z_0-ate_A),
                                                 copartisans_underestimateout=mean(Y_Z_2-Y_Z_0-ate_B),
                                                 outpartisans_underestimateout=mean(Y_Z_3-Y_Z_0-ate_C),
                                                 subset=(observed_cov5==1))

estimands_2.1.overestimate_o <- declare_inquiry(american_overestimateout=mean(Y_Z_1-Y_Z_0-ate_A),
                                                copartisans_overestimateout=mean(Y_Z_2-Y_Z_0-ate_B),
                                                outpartisans_overestimateout=mean(Y_Z_3-Y_Z_0-ate_C),
                                                subset=(observed_cov6==1))

estimands_2.2.accurate_a <- declare_inquiry(american_aa=mean(Y_Z_1-Y_Z_0-ate_A),
                                            subset=(observed_cov1==0 & observed_cov2==0))

estimands_2.2.underestimate_a <- declare_inquiry(american_ua=mean(Y_Z_1-Y_Z_0-ate_A),
                                                 subset=(observed_cov1==1))

estimands_2.2.overestimate_a <- declare_inquiry(american_oa=mean(Y_Z_1-Y_Z_0-ate_A),
                                                subset=(observed_cov2==1))

estimands_3.1 <- declare_inquiry(american=mean(Y_Z_1-Y_Z_0-ate_A),
                                 copartisans=mean(Y_Z_2-Y_Z_0-ate_B),
                                 outpartisans=mean(Y_Z_3-Y_Z_0-ate_C))

estimands_3.2 <- declare_inquiry(american=mean(Y_Z_1-Y_Z_0-ate_A))

estimands_4.dems <- declare_inquiry(american_d=mean(Y_Z_1-Y_Z_0-ate_A),
                                    copartisans_d=mean(Y_Z_2-Y_Z_0-ate_B),
                                    outpartisans_d=mean(Y_Z_3-Y_Z_0-ate_C),
                                    subset=(observed_cov1==0))

estimands_4.reps <- declare_inquiry(american_r=mean(Y_Z_1-Y_Z_0-ate_A),
                                    copartisans_r=mean(Y_Z_2-Y_Z_0-ate_B),
                                    outpartisans_r=mean(Y_Z_3-Y_Z_0-ate_C),
                                    subset=(observed_cov1==1))

estimands_5.dems <- declare_inquiry(american_d=mean(Y_Z_1-Y_Z_0-ate_A),
                                    copartisans_d=mean(Y_Z_2-Y_Z_0-ate_B),
                                    outpartisans_d=mean(Y_Z_3-Y_Z_0-ate_C),
                                    subset=(observed_cov1==0))

estimands_5.reps <- declare_inquiry(american_r=mean(Y_Z_1-Y_Z_0-ate_A),
                                    copartisans_r=mean(Y_Z_2-Y_Z_0-ate_B),
                                    outpartisans_r=mean(Y_Z_3-Y_Z_0-ate_C),
                                    subset=(observed_cov1==1))

estimands_6.nofc <- declare_inquiry(american_nofc=mean(Y_Z_1-Y_Z_0-ate_A),
                                    copartisans_nofc=mean(Y_Z_2-Y_Z_0-ate_B),
                                    outpartisans_nofc=mean(Y_Z_3-Y_Z_0-ate_C),
                                    subset=(observed_cov1==0 & observed_cov2==0))

estimands_6_w2fc <- declare_inquiry(american_w2fc=mean(Y_Z_1-Y_Z_0-ate_A),
                                    copartisans_w2fc=mean(Y_Z_2-Y_Z_0-ate_B),
                                    outpartisans_w2fc=mean(Y_Z_3-Y_Z_0-ate_C),
                                    subset=(observed_cov1==1 & observed_cov2==0))

estimands_6_w3fc <- declare_inquiry(american_w3fc=mean(Y_Z_1-Y_Z_0-ate_A),
                                    copartisans_w3fc=mean(Y_Z_2-Y_Z_0-ate_B),
                                    outpartisans_w3fc=mean(Y_Z_3-Y_Z_0-ate_C),
                                    subset=(observed_cov1==0 & observed_cov2==1))

estimands_6_bothfc <- declare_inquiry(american_bothfc=mean(Y_Z_1-Y_Z_0-ate_A),
                                      copartisans_bothfc=mean(Y_Z_2-Y_Z_0-ate_B),
                                      outpartisans_bothfc=mean(Y_Z_3-Y_Z_0-ate_C),
                                      subset=(observed_cov1==1 & observed_cov2==1))

estimands_7.nofc <- declare_inquiry(american_nofc=mean(Y_Z_1-Y_Z_0-ate_A),
                                    copartisans_nofc=mean(Y_Z_2-Y_Z_0-ate_B),
                                    outpartisans_nofc=mean(Y_Z_3-Y_Z_0-ate_C),
                                    subset=(observed_cov1==0 & observed_cov2==0))

estimands_7_w2fc <- declare_inquiry(american_w2fc=mean(Y_Z_1-Y_Z_0-ate_A),
                                    copartisans_w2fc=mean(Y_Z_2-Y_Z_0-ate_B),
                                    outpartisans_w2fc=mean(Y_Z_3-Y_Z_0-ate_C),
                                    subset=(observed_cov1==1 & observed_cov2==0))

estimands_7_w3fc <- declare_inquiry(american_w3fc=mean(Y_Z_1-Y_Z_0-ate_A),
                                    copartisans_w3fc=mean(Y_Z_2-Y_Z_0-ate_B),
                                    outpartisans_w3fc=mean(Y_Z_3-Y_Z_0-ate_C),
                                    subset=(observed_cov1==0 & observed_cov2==1))

estimands_7_bothfc <- declare_inquiry(american_bothfc=mean(Y_Z_1-Y_Z_0-ate_A),
                                      copartisans_bothfc=mean(Y_Z_2-Y_Z_0-ate_B),
                                      outpartisans_bothfc=mean(Y_Z_3-Y_Z_0-ate_C),
                                      subset=(observed_cov1==1 & observed_cov2==1))

estimands_8 <- declare_inquiry(american=mean(Y_Z_1-Y_Z_0-ate_A),
                               copartisans=mean(Y_Z_2-Y_Z_0-ate_B),
                               outpartisans=mean(Y_Z_3-Y_Z_0-ate_C))

estimands_9.dems <- declare_inquiry(american_d=mean(Y_Z_1-Y_Z_0-ate_A),
                                    copartisans_d=mean(Y_Z_2-Y_Z_0-ate_B),
                                    outpartisans_d=mean(Y_Z_3-Y_Z_0-ate_C),
                                    subset=(observed_cov1==0))

estimands_9.reps <- declare_inquiry(american_r=mean(Y_Z_1-Y_Z_0-ate_A),
                                    copartisans_r=mean(Y_Z_2-Y_Z_0-ate_B),
                                    outpartisans_r=mean(Y_Z_3-Y_Z_0-ate_C),
                                    subset=(observed_cov1==1))

## Reveal
reveal <- declare_reveal(Y,Z)

## Estimation Procedure
estimator1.1 <- declare_estimator(
  Y~as.factor(Z),
  model=lm_robust,
  term=c("as.factor(Z)1","as.factor(Z)2","as.factor(Z)3"),
  inquiry=c("american","copartisans","outpartisans")
)

estimator1.2 <- declare_estimator(
  Y~as.factor(Z),
  model=lm_robust,
  inquiry=c("american")
)

estimator2.1 <- declare_estimator(
  Y~as.factor(Z)*observed_cov1+as.factor(Z)*observed_cov2+as.factor(Z)*observed_cov3+as.factor(Z)*observed_cov4+as.factor(Z)*observed_cov5+as.factor(Z)*observed_cov6,
  model=lm_robust,
  term=c("as.factor(Z)1","as.factor(Z)1:observed_cov1","as.factor(Z)1:observed_cov2",
         "as.factor(Z)2","as.factor(Z)2:observed_cov3","as.factor(Z)2:observed_cov4",
         "as.factor(Z)3","as.factor(Z)3:observed_cov5","as.factor(Z)3:observed_cov6"),
  inquiry=c("american_accuratenat","american_underestimatenat","american_overestimatenat","copartisans_accurateco","copartisans_underestimateco","copartisans_overestimateco","outpartisans_accurateout","outpartisans_underestimateout","outpartisans_overestimateout")
)

estimator2.2 <- declare_estimator(
  Y~as.factor(Z)*observed_cov1+as.factor(Z)*observed_cov2,
  model=lm_robust,
  term=c("as.factor(Z)1","as.factor(Z)1:observed_cov1","as.factor(Z)1:observed_cov2"),
  inquiry=c("american_aa","american_ua","american_oa")
)

estimator3.1 <- declare_estimator(
  Y~as.factor(Z),
  model=lm_robust,
  term=c("as.factor(Z)1","as.factor(Z)2","as.factor(Z)3"),
  inquiry=c("american","copartisans","outpartisans")
)

estimator3.2 <- declare_estimator(
  Y~as.factor(Z),
  model=lm_robust,
  inquiry=c("american")
)

estimator4 <- declare_estimator(
  Y~as.factor(Z)*observed_cov1,
  model=lm_robust,
  term=c("as.factor(Z)1","as.factor(Z)1:observed_cov1",
         "as.factor(Z)2","as.factor(Z)2:observed_cov1",
         "as.factor(Z)3","as.factor(Z)3:observed_cov1"),
  inquiry=c("american_d","copartisans_d","outpartisans_d","american_r","copartisans_r","outpartisans_r")
)

estimator5 <- declare_estimator(
  Y~as.factor(Z)*observed_cov1,
  model=lm_robust,
  term=c("as.factor(Z)1","as.factor(Z)1:observed_cov1",
         "as.factor(Z)2","as.factor(Z)2:observed_cov1",
         "as.factor(Z)3","as.factor(Z)3:observed_cov1"),
  inquiry=c("american_d","copartisans_d","outpartisans_d","american_r","copartisans_r","outpartisans_r")
)

estimator6 <- declare_estimator(
  Y~as.factor(Z)*observed_cov1+as.factor(Z)*observed_cov2,
  model=lm_robust,
  term=c("as.factor(Z)1","as.factor(Z)1:observed_cov1","as.factor(Z)1:observed_cov2",
         "as.factor(Z)2","as.factor(Z)2:observed_cov1","as.factor(Z)2:observed_cov2",
         "as.factor(Z)3","as.factor(Z)3:observed_cov1","as.factor(Z)3:observed_cov2"),
  inquiry=c("american_nofc","american_w2fc","american_w3fc","copartisans_nofc","copartisans_w2fc","copartisans_w3fc","outpartisans_nofc","outpartisans_w2fc","outpartisans_w3fc")
)

estimator7 <- declare_estimator(
  Y~as.factor(Z)*observed_cov1+as.factor(Z)*observed_cov2,
  model=lm_robust,
  term=c("as.factor(Z)1","as.factor(Z)1:observed_cov1","as.factor(Z)1:observed_cov2",
         "as.factor(Z)2","as.factor(Z)2:observed_cov1","as.factor(Z)2:observed_cov2",
         "as.factor(Z)3","as.factor(Z)3:observed_cov1","as.factor(Z)3:observed_cov2"),
  inquiry=c("american_nofc","american_w2fc","american_w3fc","copartisans_nofc","copartisans_w2fc","copartisans_w3fc","outpartisans_nofc","outpartisans_w2fc","outpartisans_w3fc")
)

estimator8 <- declare_estimator(
  Y~as.factor(Z),
  model=lm_robust,
  term=c("as.factor(Z)1","as.factor(Z)2","as.factor(Z)3"),
  inquiry=c("american","copartisans","outpartisans")
)

estimator9 <- declare_estimator(
  Y~as.factor(Z)*observed_cov1,
  model=lm_robust,
  term=c("as.factor(Z)1","as.factor(Z)1:observed_cov1",
         "as.factor(Z)2","as.factor(Z)2:observed_cov1",
         "as.factor(Z)3","as.factor(Z)3:observed_cov1"),
  inquiry=c("american_d","copartisans_d","outpartisans_d","american_r","copartisans_r","outpartisans_r")
)

set.seed(1000)
design1.1 <- population1.1+potential_1.1+estimands_1.1+randomization_partisans+reveal+estimator1.1
set.seed(1000)
design1.2 <- population1.2+potential_1.2+estimands_1.2+randomization_full+reveal+estimator1.2
set.seed(1000)
design2.1 <- population2.1+potential_2.1+estimands_2.1.accurate_a+randomization_partisans+reveal+estimator2.1
set.seed(1000)
design2.2 <- population2.2+potential_2.2+estimands_2.2.accurate_a+randomization_full+reveal+estimator2.2
set.seed(1000)
design3.1 <- population3.1+potential_3.1+estimands_3.1+randomization_partisans+reveal+estimator3.1
set.seed(1000)
design3.2 <- population3.2+potential_3.2+estimands_3.2+randomization_full+reveal+estimator3.2
set.seed(1000)
design4 <- population4+potential_4+estimands_4.reps+randomization_partisans+reveal+estimator4
set.seed(1000)
design5 <- population5+potential_5+estimands_5.reps+randomization_partisans+reveal+estimator5
set.seed(1000)
design6 <- population6+potential_6+estimands_6_w2fc+randomization_partisans+reveal+estimator6
set.seed(1000)
design7 <- population7+potential_7+estimands_7_w2fc+randomization_partisans+reveal+estimator7
set.seed(1000)
design8 <- population8+potential_8+estimands_8+randomization_partisans+reveal+estimator8
set.seed(1000)
design9 <- population9+potential_9+estimands_9.reps+randomization_partisans+reveal+estimator9

## Simulations
set.seed(1000)
d1.1 <- diagnose_design(design1.1,diagnosands=declare_diagnosands(power=mean(p.value<.05)),sims=500)
set.seed(1000)
d1.2 <- diagnose_design(design1.2,diagnosands=declare_diagnosands(power=mean(p.value<.05)),sims=500)
set.seed(1000)
d2.1 <- diagnose_design(design2.1,diagnosands=declare_diagnosands(power=mean(p.value<.05)),sims=500)
set.seed(1000)
d2.2 <- diagnose_design(design2.2,diagnosands=declare_diagnosands(power=mean(p.value<.05)),sims=500)
set.seed(1000)
d3.1 <- diagnose_design(design3.1,diagnosands=declare_diagnosands(power=mean(p.value<.05)),sims=500)
set.seed(1000)
d3.2 <- diagnose_design(design3.2,diagnosands=declare_diagnosands(power=mean(p.value<.05)),sims=500)
set.seed(1000)
d4 <- diagnose_design(design4,diagnosands=declare_diagnosands(power=mean(p.value<.05)),sims=500)
set.seed(1000)
d5 <- diagnose_design(design5,diagnosands=declare_diagnosands(power=mean(p.value<.05)),sims=500)
set.seed(1000)
d6 <- diagnose_design(design6,diagnosands=declare_diagnosands(power=mean(p.value<.05)),sims=500)
set.seed(1000)
d7 <- diagnose_design(design7,diagnosands=declare_diagnosands(power=mean(p.value<.05)),sims=500)
set.seed(1000)
d8 <- diagnose_design(design8,diagnosands=declare_diagnosands(power=mean(p.value<.05)),sims=500)
set.seed(1000)
d9 <- diagnose_design(design9,diagnosands=declare_diagnosands(power=mean(p.value<.05)),sims=500)

# Detect 80% power on interactions, Main Text
## Footnote 7, American x underestimated
set.seed(1000)
redesign2.1.a1 <- redesign(design2.1,ate_b_2.1=seq(0,1,by=0.01))
set.seed(1000)
redesign2.1.a1.diag <- diagnose_design(redesign2.1.a1,diagnosands = declare_diagnosands(power = mean(p.value<.05)),sims=500) ## 0.32 to reach 80% power for americans/underestimators
subset(redesign2.1.a1.diag$diagnosands_df,term=="as.factor(Z)2")
## Footnote 8, Co-Partisan x underestimated
set.seed(1000)
redesign2.1.c1 <- redesign(design2.1,ate_a_cov2_2.1=seq(0,1,by=0.01)) ## power=0.50
set.seed(1000)
redesign2.1.c1.diag <- diagnose_design(redesign2.1.c1,diagnosands = declare_diagnosands(power = mean(p.value<.05)),sims=500) ## 0.13 to reach 80% power
subset(redesign2.1.c1.diag$diagnosands_df,term=="as.factor(Z)1:observed_cov2")
## Footnote 9, Out-Partisan x underestimated
set.seed(1000)
redesign2.1.o1 <- redesign(design2.1,ate_c_cov5_2.1=seq(0,1,by=0.01)) ## power=0.32
set.seed(1000)
redesign2.1.o1.diag <- diagnose_design(redesign2.1.o1,diagnosands = declare_diagnosands(power = mean(p.value<.05)),sims=500) ## 0.13 to reach 80% power
subset(redesign2.1.o1.diag$diagnosands_df,term=="as.factor(Z)3:observed_cov5")
## Footnote 10, American x Republican (Wearing)
set.seed(1000)
redesign4.a1 <- redesign(design4,ate_a_cov1_4=seq(0,0.6,by=0.01))
set.seed(1000)
redesign4.a1.diag <- diagnose_design(redesign4.a1,diagnosands=declare_diagnosands(power=mean(p.value<.05)),sims=500)
subset(redesign4.a1.diag$diagnosands_df,term=="as.factor(Z)1:observed_cov1")
## Footnote 10, Co-Partisan x Republican (Wearing)
set.seed(1000)
redesign4.c1 <- redesign(design4,ate_b_cov1_4=seq(0,0.6,by=0.01))
set.seed(1000)
redesign4.c1.diag <- diagnose_design(redesign4.c1,diagnosands=declare_diagnosands(power=mean(p.value<.05)),sims=500)
subset(redesign4.c1.diag$diagnosands_df,term=="as.factor(Z)2:observed_cov1")
## Footnote 10, Out-Partisan x Republian (Wearing)
set.seed(1000)
redesign4.o1 <- redesign(design4,ate_c_cov1_4=seq(0,0.6,by=0.01))
set.seed(1000)
redesign4.o1.diag <- diagnose_design(redesign4.o1,diagnosands=declare_diagnosands(power=mean(p.value<.05)),sims=500)
subset(redesign4.o1.diag$diagnosands_df,term=="as.factor(Z)3:observed_cov1")
## Footnote 12, American x Republican (Effectiveness)
set.seed(1000)
redesign5.a1 <- redesign(design5,ate_a_cov1_5=seq(0,0.6,by=0.01))
set.seed(1000)
redesign5.a1.diag <- diagnose_design(redesign5.a1,diagnosands=declare_diagnosands(power=mean(p.value<.05)),sims=500)
subset(redesign5.a1.diag$diagnosands_df,term=="as.factor(Z)1:observed_cov1")
## Footnote 12, Co-Partisan x Republican (Effectiveness)
set.seed(1000)
redesign5.c1 <- redesign(design5,ate_b_cov1_5=seq(0,0.6,by=0.01))
set.seed(1000)
redesign5.c1.diag <- diagnose_design(redesign5.c1,diagnosands=declare_diagnosands(power=mean(p.value<.05)),sims=500)
subset(redesign5.c1.diag$diagnosands_df,term=="as.factor(Z)2:observed_cov1")
## Footnote 12, Out-Partisan x Republican (Effectiveness)
set.seed(1000)
redesign5.o1 <- redesign(design5,ate_c_cov1_5=seq(0,0.6,by=0.01))
set.seed(1000)
redesign5.o1.diag <- diagnose_design(redesign5.o1,diagnosands=declare_diagnosands(power=mean(p.value<.05)),sims=500)
subset(redesign5.o1.diag$diagnosands_df,term=="as.factor(Z)3:observed_cov1")

# Detect 80% Power on Interactions, Appendix
## Footnote 1, American x W2
set.seed(1000)
redesign6.c <- redesign(design6,ate_b_6=seq(0,0.6,by=0.01))
set.seed(1000)
redesign6.c.diag <- diagnose_design(redesign6.c,diagnosands=declare_diagnosands(power=mean(p.value<.05)),sims=500)
subset(redesign6.c.diag$diagnosands_df,term=="as.factor(Z)2")
## Footnote 1, American x W3
set.seed(1000)
redesign6.o <- redesign(design6,ate_c_6=seq(0,0.6,by=0.01))
set.seed(1000)
redesign6.o.diag <- diagnose_design(redesign6.o,diagnosands=declare_diagnosands(power=mean(p.value<.05)),sims=500)
subset(redesign6.o.diag$diagnosands_df,term=="as.factor(Z)3")
## Footnote 2, Co-Partisan x W3
set.seed(1000)
redesign6.o1 <- redesign(design6,ate_c_cov1_6=seq(0,0.6,by=0.01))
set.seed(1000)
redesign6.o1.diag <- diagnose_design(redesign6.o1,diagnosands=declare_diagnosands(power=mean(p.value<.05)),sims=500)
subset(redesign6.o1.diag$diagnosands_df,term=="as.factor(Z)3:observed_cov1")